home *** CD-ROM | disk | FTP | other *** search
/ Our Solar System / Our Solar System.iso / shuttle / seesat3 / astro.c next >
Encoding:
C/C++ Source or Header  |  1990-12-24  |  17.0 KB  |  630 lines
  1. /*
  2. ASTRO.C
  3. by Paul S. Hirose, 1990 Nov 12
  4. Coordinate transformation and astronomical functions for SEESAT satellite
  5. tracking program.
  6.  
  7. This file is in the public domain.
  8.  
  9. 200 - 299
  10. */
  11.  
  12. #include "B:SEESAT.H"    /* global header */
  13.  
  14. /* library functions:
  15.     asin atan atan2 cos fabs printf sin sqrt tan
  16. are used in this file. */
  17.  
  18.  
  19. extern void printf();
  20. extern double asin(), atan(), atan2(), cos(), fabs(),
  21.     sin(), sqrt(), tan();
  22.  
  23. extern double dint();
  24.  
  25. /* dint(3.6) == 3.0, dint(-3.6) == -3.0.
  26. If your library has floor() and ceil() instead of dint(), comment out the
  27. preceding dint() declaration and un-comment the following function: */
  28.  
  29. /*
  30. static double dint(x)
  31. double x;
  32. {
  33.     extern double floor(), ceil();
  34.  
  35.     if (x >= 0.)
  36.         return floor(x);
  37.     else
  38.         return ceil(x);
  39. }
  40. */
  41.  
  42. /* Default observer's location.  All four quantities may be expressed to any
  43. desired number of decimal places. */
  44. #define D_ALT "2400"        /* FEET above sea level */
  45. #define D_LAT "34.9147"        /* north latitude */
  46. #define D_LON "-117.9333"    /* east longitude */
  47. #define D_TZ "8"        /* UTC - local time, hours */
  48.  
  49. #define FLAT 3.35278e-3        /* flattening of earth.
  50.                 WGS-72 value = 3.35278e-3 */
  51.  
  52. /* REFEP is the default epoch to which R.A./dec. will be precessed.
  53. E.g. (per Meeus p. 101):
  54. 3477629251. = epoch 1900.0, 3503926689. = 1950.0, 3530224800. = 2000.0.
  55. EPMSG is part of the message (printed at startup) that states the epoch for
  56. R.A./dec.  If you've enabled precession, be sure REFEP & EPMSG agree. */
  57.  
  58. #if ENPRE        /* precession enabled */
  59. #define REFEP 3530224800.
  60. #define EPMSG "2000.0"
  61.  
  62. #else            /* don't precess */
  63. #define EPMSG "= epoch of elements"
  64. #endif
  65.  
  66.  
  67. /* If GLON is nonzero, the displayed satellite longitude is its longitude
  68. relative to Greenwich.  If zero, longitude is relative to observer's
  69. meridian.  GLON only affects satellite longitude readout. */
  70.  
  71. #define GLON 0
  72.  
  73. #if GLON
  74. #define LMSG "Greenwich"
  75. #else
  76. #define LMSG "observer's meridian"
  77. #endif
  78.  
  79.  
  80. /* 0 for normal operation.  Non-zero will make static functions and data
  81. extern */
  82. #define TEST 0
  83.  
  84.  
  85. /* ALL should normally be 0, which will suppress printing of predictions when
  86. satellite is below horizon.  To see all predictions regardless of whether or
  87. not satellite is visible, set ALL to 1. */
  88.  
  89. #define ALL 0
  90.  
  91. /*######################## LOCAL FUNCTIONS ########################*/
  92.  
  93. #if TEST
  94. #define STATIC        /* precedes static func definitions */
  95. #define SC extern    /* precedes static func declarations */
  96.  
  97. #else
  98. #define STATIC static
  99. #define SC static
  100. #endif
  101.  
  102. SC void eqhor();        /* convert rectangular to az/el, print */
  103. SC void eceq();            /* ecliptical to equatorial coordinates */
  104. SC void point();        /* align +z axis with a vector */
  105. SC void recsph();        /* rectangular to spherical coordinates */
  106. SC void sun();            /* xyz of sun (interpol.) */
  107. SC double sunlon();        /* mean longitude of sun */
  108. SC void sunref();        /* xyz of sun (non-interpol.) */
  109. SC void yrot();            /* y-rotate vector */
  110. SC void zrot();            /* z-rotate vector */
  111.  
  112. /*########################### LOCAL DATA ###########################*/
  113.  
  114. #if ENPRE
  115.  
  116. /* dircos[] is a rotation matrix.  Initialized by inpre() & used to precess
  117. satellite Right Ascension & declination.  terep is the epoch to which the
  118. coordinates will be precessed. */
  119.  
  120. STATIC struct {
  121.     double ix, iy, iz, jx, jy, jz, kx, ky, kz;
  122. } dircos;
  123. STATIC double terep = REFEP;
  124.  
  125. #endif
  126.  
  127. /* data for observer's location */
  128. STATIC struct {
  129.     double
  130.     h,        /* altitude above sea level */
  131.     lambda,        /* longitude */
  132.     sinlat, coslat,    /* latitude sin & cos */
  133.     zg,        /* north distance, observer to equatorial plane */
  134.     xc;        /* distance, observer to polar axis */
  135. } obs;
  136.  
  137. /* Sin & cos of epsilon, the obliquity of ecliptic.  Since epsilon changes
  138. but .013 deg/century, it's sufficient to use a fixed value, that of 1975
  139. (23.442 deg). */
  140. STATIC double sineps = .39783;
  141. STATIC double coseps = .91746;
  142.  
  143. /*############################## CODE ##############################*/
  144.  
  145. STATIC void
  146. eqhor(equa, t)
  147. struct vector *equa;    /* vector to object */
  148. double t;        /* epoch */
  149. /* Converts vector *equa from a rectangular (Aries, equator, north) system to
  150. a polar (azimuth, elevation, range) system at epoch t.  Prints the azimuth &
  151. elevation rounded to nearest deg.  No correction for refraction or geocentric
  152. parallax. */
  153. {
  154.     double    lhaa;    /* local hour angle of Aries */
  155.  
  156.     /* rotate to a south, east, zenith system */
  157.     lhaa = thetag(t) + obs.lambda;
  158.     zrot(equa, sin(lhaa), cos(lhaa));
  159.     yrot(equa, obs.coslat, obs.sinlat);
  160.  
  161.     equa->x = -equa->x;    /* makes az run correct direction */
  162.     recsph(equa);        /* convert to spherical coordinates */
  163.  
  164.     printf("az=%d ", (int) (equa->x * ra2de + .5));
  165.     printf("el=%d\n", (int) (equa->y * ra2de +
  166.       (equa->y >= 0. ? .5 : -.5)));
  167. }
  168.  
  169.  
  170. void
  171. dusk()
  172. /* Print azimuth and elevation of the sun.  The next argument(s) on
  173. the command line are taken as the epoch. */
  174. {
  175.     double t;
  176.     struct vector csun;    /* equatorial coordinates of sun */
  177.  
  178.     t = tokmin();        /* get epoch from cmd line */
  179.     sun(&csun, t);        /* get sun's position */
  180.     eqhor(&csun, t);    /* convert to az/el, print */
  181. }
  182.  
  183.  
  184. STATIC void
  185. eceq(eclip)
  186. struct vector *eclip;    /* ecliptical coordinates */
  187. /* Coverts *eclip to equatorial coordinates, i.e., x-rotates by the negative
  188. obliquity of the ecliptic. */
  189. {
  190.     double tempy;
  191.  
  192.     tempy    = eclip->y * coseps - eclip->z * sineps;
  193.     eclip->z = eclip->y * sineps + eclip->z * coseps;
  194.     eclip->y = tempy;
  195. }
  196.  
  197.  
  198. double
  199. fmod2p(x)
  200. double x;
  201. /* Reduces x to range 0 - 2pi */
  202. {
  203.     x /= twopi;
  204.     x = (x - dint(x)) * twopi;
  205.     if (x < 0.)
  206.         return x + twopi;
  207.     else
  208.         return x;
  209. }
  210.  
  211.  
  212. #if ENPRE    /* precession code enabled */
  213.  
  214. void
  215. inpre()
  216. /* Initialize rotation matrix dircos to precess from "epoch" (a global
  217. variable representing epoch of the satellite elements) to terep.  terep is an
  218. external variable in this file.  Assumes that the effect of precession is a
  219. uniform 50".4 increase of longitude per year. */
  220. {
  221.     double a, sina, cosa;
  222.     struct vector *vecp;
  223.     int i;
  224.  
  225.     a = (terep - epoch) * -4.65e-10;    /* luni-solar precession */
  226.     sina = sin(a);
  227.     cosa = cos(a);
  228.  
  229.     /* Ecliptical components of (I, J, K) unit vectors */
  230.     dircos.ix = 1.;        dircos.iy = 0.;        dircos.iz = 0.;
  231.     dircos.jx = 0.;        dircos.jy = coseps;    dircos.jz = -sineps;
  232.     dircos.kx = 0.;        dircos.ky = sineps;    dircos.kz = coseps;
  233.  
  234.     /* Rotate each vector in longitude by the amount of luni-solar
  235.     precession.  Then convert from ecliptical to equatorial coordinates.
  236.     We end up with the old unit vectors expressed in terms of the new
  237.     (i.e., precessed) equatorial coordinate system. */
  238.  
  239.     for (vecp = (struct vector *) &dircos, i = 3; i--; ++vecp) {
  240.         zrot(vecp, sina, cosa);
  241.         eceq(vecp);
  242. }    }
  243.  
  244. #endif
  245.  
  246.  
  247. void
  248. moon()
  249. /* Prints the azimuth, elevation, and % of illumination of the moon.  Epoch
  250. is taken from the next arguments on the command line.  Position is good
  251. within a couple degrees or so. */
  252. {
  253.     struct vector luna, sol, terra;
  254.     double t, t1900, lp, mp, f, lamb, beta, cosb;
  255.  
  256.     t = tokmin();            /* obtain epoch from cmd line */
  257.  
  258.     /* formulas from Meeus */
  259.     t1900 = t - 3477628800.;        /* 1900 Jan 0 12h UT */
  260.     lp = fmod2p(4.72 + 1.5970243e-4 * t1900);    /* longitude */
  261.     mp = fmod2p(5.168 + 1.5835217e-4 * t1900);    /* mean anomaly */
  262.     f = fmod2p(.1964 + 1.6034425e-4 * t1900);    /* dist fm asc node */
  263.     lamb = lp + .1097 * sin(mp);            /* longitude */
  264.     beta = .0895 * sin(f);                /* latitude */
  265.  
  266.     /* get ecliptical components of unit vector to moon */
  267.     luna.z = sin(beta);
  268.     cosb = cos(beta);
  269.     luna.y = sin(lamb) * cosb;
  270.     luna.x = cos(lamb) * cosb;
  271.  
  272.     eceq(&luna);        /* ecliptical to equatorial */
  273.  
  274.     /* terra = unit vector moon to earth */
  275.     terra.x = -luna.x;
  276.     terra.y = -luna.y;
  277.     terra.z = -luna.z;
  278.  
  279.     eqhor(&luna, t);    /* print az el */
  280.  
  281.     /* Get moon->sun unit vector.  sun(), which yields the earth->sun
  282.     vector, comes close enough.   Rotate axes of terra to point +z axis
  283.     at sun.  This will make terra.z equal to the cos of the moon's phase
  284.     angle.  Plug that into Meeus' illumination formula, print. */
  285.  
  286.     sun(&sol, t);
  287.     point(&terra, &sol);
  288.  
  289.     printf("%d%% illuminated\n", (int) (50. * (1. + terra.z)));
  290. }
  291.  
  292.  
  293. void
  294. parall()
  295. /* Prints the parallactic angle:  the angle (at the satellite) between
  296. celestial north and observer's zenith.  The next argument(s) on the command
  297. line are taken as the desired epoch. */
  298. {
  299.     struct vector zenith, topsat;
  300.     double lhaa, t, coslha, sinlha;
  301.  
  302.     t = tokmin();
  303.     lhaa = thetag(t) + obs.lambda;    /* local hr angle of Aries */
  304.     sinlha = sin(lhaa);
  305.     coslha = cos(lhaa);
  306.  
  307.     /* Get equatorial coordinates of satellite at time t. */
  308.     sgp4(t - epoch - toffs);
  309.  
  310.     /* get topocentric sat coordinates */
  311.     topsat.x = sat.x - obs.xc * coslha;
  312.     topsat.y = sat.y - obs.xc * sinlha;
  313.     topsat.z = sat.z + obs.zg;
  314.  
  315.     /* Load struct "zenith" with a unit vector directed from the observer
  316.     to the zenith, expressed in equatorial coordinates. */
  317.  
  318.     zenith.z = obs.sinlat;
  319.     zenith.x = coslha * obs.coslat;
  320.     zenith.y = sinlha * obs.coslat;
  321.  
  322.     /* z-rotate the axes of "zenith" so satellite is in the x-z plane of
  323.     "zenith" (+x half plane).  Then y-rotate axes to point +z axis at
  324.     satellite. */
  325.  
  326.     point(&zenith, &topsat);
  327.  
  328.     /* The +z axis passes through the satellite.  The north pole
  329.     has a y-coordinate of 0 and a negative x-coordinate.  The
  330.     origin of the coordinate system is still at the observer. */
  331.  
  332.     printf("parallactic angle = %d\n",
  333.       (int) (fmod2p(atan2(zenith.y, -zenith.x)) * ra2de));
  334. }
  335.  
  336.  
  337. STATIC void
  338. point(v, d)
  339. struct vector *v, *d;
  340. /* rotates the axes of vector v to align the +z axis with vector d. */
  341. {
  342.     double fz, fz2, r;
  343.  
  344.     fz2 = d->x * d->x + d->y * d->y;
  345.     fz = sqrt(fz2);
  346.     zrot(v, d->y / fz, d->x / fz);
  347.     r = sqrt(fz2 + d->z * d->z);
  348.     yrot(v, fz / r, d->z / r);
  349. }
  350.  
  351.  
  352. STATIC void
  353. recsph(v)
  354. struct vector *v;
  355. /* Converts vector v from rectangular to spherical coordinates.  On exit,
  356. v->x = angle in the x-y plane, v->y = angle with respect to the x-y plane,
  357. and v->z = radius.  Angle x is 0 to 2pi, y is -pi/2 to pi/2. */
  358. {
  359.     double temp;
  360.  
  361.     temp = sqrt(v->x * v->x + v->y * v->y + v->z * v->z);
  362.     v->x = atan2(v->y, v->x);
  363.     if (v->x < 0.)
  364.         v->x += twopi;
  365.     v->y = asin(v->z / temp);
  366.     v->z = temp;
  367. }
  368.  
  369.  
  370. #if ENPRE    /* precession code enabled */
  371.  
  372. void
  373. setep()
  374. /* Sets rotation matrix "dircos" to precess from "epoch" (a global variable
  375. representing epoch of the satellite orbital elements) to the date/time given
  376. on the command line. */
  377. {
  378.     terep = tokmin();    /* get epoch from cmd line */
  379.     inpre();    /* reinitialize rotation matrix */
  380. }
  381.  
  382. #endif
  383.  
  384.  
  385. STATIC void
  386. sun(pos, ep)
  387. struct vector *pos;
  388. double ep;
  389. {
  390. /* Equatorial xyz components of a unit vector to sun at epoch "ep" are
  391. returned in struct "pos".  Accuracy about .1 deg. */
  392.  
  393.     static double t0,    /* start epoch of a 6-day interval */
  394.         dx, dy, dz;    /* change in x, y, z during 6 days */
  395.     static struct vector sun0;    /* sun position at t0 */
  396.  
  397.     double t;    /* e.g. .5 if ep = t0 + 3 days */
  398.  
  399.     if ((t = (ep - t0) / 8640.) > 1.0  ||  t < 0.) {
  400.     /* Position was requested for a time outside the current 6-day
  401.     interval.  So set up a new 6-day interval centered at "ep". */
  402.         struct vector sun1;
  403.  
  404.         t0 = ep - 4320.;    /* 3 days before "ep" */
  405.         sunref(&sun0, t0);
  406.         sunref(&sun1, ep + 4320.);    /* 3 days after "ep" */
  407.         dx = sun1.x - sun0.x;
  408.         dy = sun1.y - sun0.y;
  409.         dz = sun1.z - sun0.z;
  410.         t = .5;
  411.     }
  412.     /* get position by linear interpolation in the 6-day period */
  413.     pos->x = sun0.x + t * dx;
  414.     pos->y = sun0.y + t * dy;
  415.     pos->z = sun0.z + t * dz;
  416. }
  417.  
  418.  
  419. STATIC double
  420. sunlon(epoch)
  421. double epoch;
  422. /* Return geometric mean longitude of sun at "epoch", accurate to about .1
  423. deg.  Formulas adapted from Meeus.  The returned value may be off by a small
  424. multiple of 2pi.  That should cause no problem since in this program the
  425. longitude is only used by sin() & cos(). */
  426. {
  427.     static double e = .016709;    /* eccentricity year 2000 */
  428.  
  429.     double t, l, m, enew, eold, v;
  430.  
  431.     t = epoch - 3477628800.;    /* since 1900 Jan 0 12h */
  432.  
  433.     /* geometric mean long. */
  434.     l = fmod2p(4.881628 + t * 1.1946383e-5);
  435.  
  436.     /* mean anomaly */
  437.     m = fmod2p(6.256584 + t * 1.1945812e-5);
  438.  
  439.     /* solve Kepler's equation to .1 deg. */
  440.     enew = m;
  441.     do {
  442.         eold = enew;
  443.         enew = m + e * sin(eold);
  444.     } while (fabs(enew - eold) >= 1.7e-3);
  445.  
  446.     /* tan is asymptotic at pi/2, so if eccentric anomaly is within .1
  447.     deg of pi radians, we'll just make true anomaly = pi. */
  448.     if (fabs(enew - pi) <= 1.7e-3)
  449.         v = pi;
  450.     else
  451.         v = 2. * atan(sqrt((1. + e) / (1. - e)) * tan(enew / 2.));
  452.  
  453.     return l + v - m;
  454. }
  455.  
  456.  
  457. STATIC void
  458. sunref(pos, ep)
  459. struct vector *pos;
  460. double ep;
  461. /* Put the mean equatorial xyz components of a unit vector to sun at
  462. epoch "ep" into struct "pos".  Good to about .1 deg. */
  463. {
  464.     double lon, sinlon;
  465.  
  466.     lon = sunlon(ep);
  467.     sinlon = sin(lon);
  468.  
  469.     pos->x = cos(lon);
  470.     pos->y = sinlon * coseps;
  471.     pos->z = sinlon * sineps;
  472. }
  473.  
  474.  
  475. double
  476. thetag(ep)
  477. double ep;    /* epoch */
  478. /* Returns Greenwich hour angle of the mean equinox at ep.  If ep is UT1,
  479. returned value is within .001 sec of time compared to the formula in the '85
  480. Astronomical Almanac, for 1990 through 2010.  As a side effect, this function
  481. sets ds50 to days since 1950 Jan 0 0h UT. */
  482. {
  483.     double theta;
  484.  
  485.     ds50 = (ep - 3503925360.) / xmnpda;
  486.     theta = .27524987 + 1.00273790935 * ds50;    /* revolutions */
  487.  
  488.     return (theta - dint(theta)) * twopi;
  489. }
  490.  
  491.  
  492. double
  493. topos()
  494. {
  495.     double lat, c, s;
  496.     static double f = FLAT;        /* flattening of earth */
  497.  
  498.     printf("RA & dec epoch %s,\n", EPMSG);
  499.     printf("satellite longitude relative to %s\n\n", LMSG);
  500.  
  501.     printf("observer's location:\n");
  502.     printf("feet above sea level ");
  503.     obs.h = din(D_ALT) * 4.778826e-8;    /* ft->earth radii */
  504.     printf("north latitude ");
  505.     lat = din(D_LAT) * de2ra;
  506.     printf("east longitude ");
  507.     obs.lambda = din(D_LON) * de2ra;
  508.  
  509.     obs.sinlat = sin(lat);
  510.     obs.coslat = cos(lat);
  511.  
  512.     /* observer's position with respect to geocenter (formulas
  513.     from Baker & Makemson) */
  514.     c = 1. / sqrt(1. - (f + f - f * f) * obs.sinlat * obs.sinlat);
  515.     s = c * (1. - f) * (1. - f);
  516.     obs.xc = (c + obs.h) * obs.coslat;
  517.     obs.zg = -(s + obs.h) * obs.sinlat;
  518.  
  519.     printf("UTC - local time, hours ");
  520.     return din(D_TZ) * 60.0;
  521. }
  522.  
  523.  
  524. int
  525. xyztop(t)
  526. double t;    /* epoch */
  527. /* If satellite is below horizon and ALL is #defined 0, xyztop() returns 0. 
  528. If satellite is above horizon or ALL is #defined nonzero, loads globals
  529. elsusa, azel, radec and latlon with appropriate values and returns 1. */
  530. {
  531.     struct vector suneq;    /* equatorial sun position */
  532.  
  533.     double temp, tempx, tempy, lhaa, sinlha, coslha;
  534.  
  535.     lhaa = thetag(t) + obs.lambda;        /* local hr angle Aries */
  536.     sinlha = sin(lhaa);
  537.     coslha = cos(lhaa);
  538.  
  539.     /* load azel & radec with sat coordinates translated to observer */
  540.     radec.x = azel.x = sat.x - obs.xc * coslha;
  541.     radec.y = azel.y = sat.y - obs.xc * sinlha;
  542.     radec.z = azel.z = sat.z + obs.zg;
  543.  
  544.     /* Rotate azel into a south/east/zenith system */
  545.     zrot(&azel, sinlha, coslha);
  546.     yrot(&azel, obs.coslat, obs.sinlat);
  547. #if ALL == 0
  548.     if (azel.z < 0.)
  549.         return 0;    /* below horizon */
  550. #endif
  551.     azel.x = -azel.x;    /* makes az go correct direction */
  552.     FTEST((200, 3, 4, &azel.x, &azel.y, &azel.z));
  553.     recsph(&azel);        /* convert to spherical */
  554.  
  555. #if ENPRE    /* precess RA & dec */
  556.     tempx   = radec.x * dircos.ix + radec.y * dircos.jx +
  557.       radec.z * dircos.kx;
  558.     tempy   = radec.x * dircos.iy + radec.y * dircos.jy +
  559.       radec.z * dircos.ky;
  560.     radec.z = radec.x * dircos.iz + radec.y * dircos.jz +
  561.       radec.z * dircos.kz;
  562.  
  563.     radec.x = tempx;
  564.     radec.y = tempy;
  565. #endif
  566.     FTEST((201, 3, 4, &radec.x, &radec.y, &radec.z));
  567.     recsph(&radec);        /* convert to spherical */
  568.  
  569.     /* latitude & longitude */
  570.     latlon.x = sat.x;
  571.     latlon.y = sat.y;
  572.     latlon.z = sat.z;
  573.     recsph(&latlon);    /* convert to spherical coordinates */
  574.  
  575. #if GLON    /* figure longitude relative to Greenwich */
  576.     latlon.x = fmod2p(latlon.x - lhaa + obs.lambda);
  577. #else        /* figure longitude relative to observer */
  578.     latlon.x = fmod2p(latlon.x - lhaa);
  579. #endif
  580.     if (latlon.x > pi)
  581.         latlon.x -= twopi;
  582.  
  583.     latlon.z -= mean_r;    /* mean_r = mean radius of earth */
  584.  
  585.     /* Sun elevation at satellite.  Get sun position, rotate the
  586.     axes so the satellite is on the positive z-axis.  Sun elevation =
  587.     arc sin z.  Add dip of horizon due to height of satellite. */
  588.  
  589.     sun(&suneq, t);
  590.     point(&suneq, &sat);
  591.  
  592.     temp = (asin(suneq.z) + atan(sqrt((2. + latlon.z) * latlon.z)))
  593.       * ra2de;
  594.     FTEST((202, 1, 1, &temp));
  595.     elsusa = temp + ((temp >= 0.) ? .5 : -.5);
  596.  
  597.     return 1;    /* signifies sat elev >= 0 */
  598. }
  599.  
  600.  
  601. STATIC void
  602. yrot(vec, sint, cost)
  603. struct vector *vec;
  604. double sint, cost;
  605. /* y-rotates coordinate axes for "vec" by angle t, where "sint"
  606. and "cost" are the sine and cosine of t. */
  607. {
  608.     double tempz;
  609.     tempz  = vec->x * sint + vec->z * cost;
  610.     vec->x = vec->x * cost - vec->z * sint;
  611.     vec->z = tempz;
  612. }
  613.  
  614.  
  615. STATIC void
  616. zrot(vec, sint, cost)
  617. struct vector *vec;
  618. double sint, cost;
  619. /* similar to yrot(), except rotate about z axis */
  620. {
  621.     double tempx;
  622.  
  623.     tempx  = vec->y * sint + vec->x * cost;
  624.     vec->y = vec->y * cost - vec->x * sint;
  625.     vec->x = tempx;
  626. }
  627. h */
  628.  
  629.     /* Sun elevation at satellite.  Get sun position, rotate the
  630.     axes so the s